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Abstract. We present STARDUST, a new self-consistent 
modelling of the spectral energy distributions (SEDs) of galax- 
ies from far-UV to radio wavelengths. In order to derive the 
SEDs in this broad spectral range, we first couple spectropho- 
tometric and (closed-box) chemical evolutions to account for 
metallicity effects on the spectra of synthetic stellar popula- 
tions. We briefly compare the UV/visible/near-IR colours and 
magnitudes predicted by our code with those of other codes 
available in the literature and we find an overall agreement, 
in spite of differences in the stellar data. We then use a phe- 
nomenological fit for the metal-dependent extinction curve and 
a simple geometric distribution of the dust to compute the opti- 
cal depth of galaxies and the corresponding obscuration curve. 
This enables us to calculate the fraction of stellar light re- 
processed in the infrared range. In a final step, we define a 
dust model with various components and we fix the weights 
of these components in order to reproduce the IRAS correla- 
tion of IR colours with total IR luminosities. This allows us 
to compute far-IR SEDs that phenomenologically mimic ob- 
served trends. We are able to predict the spectral evolution of 
galaxies in a broad wavelength range, and we can reproduce 
the observed SEDs of local spirals, starbursts, luminous in- 
frared galaxies (LIRGs) and ultra luminous infrared galaxies 
(ULIRGs). This modelling is so far kept as simple as possi- 
ble and depends on a small number of free parameters, namely 
the initial mass function (IMF), star formation rate (SFR) time 
scale, gas density, and galaxy age, as well as on more refined 
assumptions on dust properties and the presence (or absence) 
of gas inflows/outflows. However, these SEDs will be subse- 
quently implemented in a semi-analytic approach of galaxy 
formation, where most of the free parameters can be consis- 
tently computed from more general assumptions for the physi- 
cal processes ruling galaxy formation and evolution. 



1. Introduction 

Our knowledge of galaxies is essentially based upon the light 
they emit, so that anyone interested in studying such objects has 
to address the obvious — but nevertheless complex — issue of 
their spectrophotometric evolution. This becomes even more 
difficult if one is interested in handling their luminosity budget 
properly. Part of the light emitted by stars is reprocessed in the 
mid/far-IR window by the dust these very same stars produce, 
so that our knowledge of the star formation history of an indi- 
vidual galaxy (and a fortiori the overall star formation history 
of the universe) crucially depends on our ability to estimate the 
shape of its spectrum from the far-UV to the submm wave- 
length range. The purpose of this paper is to compute SEDs in 
the broadest wavelength range, that can be used to analyze the 
increasing amount of multiwavelength observations targetting 
high-z galaxies. 

As a matter of fact, several pieces of evidence are now con- 
verging to give a novel view of forming galaxies' luminosity 
budget. This can be summarized as follows. 

The discovery of the cosmic infrared background (CIRB) 
by the COBE satellite, is the first evidence that high-z galax- 
ies strongly emit in the IR/submm range, and that dust played 
an important role in the luminosity bud get of these galaxies 
(Puge t et al. ( |1996| ); G uiderd oni et al. ( |1997|); Fix sen et al. 
( |1998| ), Schlegel et al. ( |1998| ), Hauser et al. ( |1998| )). Though 
the exact level of the CIRB is still a matter of debate (Lagache 
et al., 1999), it appears to be 5-10 times higher than the no- 
evolution prediction based on the IRAS local luminosity func- 
tion, and 2 times more intense than the cosmic optical back- 
ground dPozzetti et al., 1998| ). 

At 850 /im, the background has been broken into its com- 
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ponent s by S CUBA at the J CMT ( Small et al. ( |1997[), H ughes 
( [199^ ), Barger et al. ( [l99^ ), Eales et al. ( |1998| )). The 



galaxy evo- same work has been achieved respectively by ISOPHOT at 



175 /im (Kawara et al. (1998), Puget et al. (1999), Dole et al. 
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( |1999|)), a n d ISO CAM at 15 /im (Oliver et al. ( |1997[ ), Elbaz et 
al. ( |1998[ ),( |l999| )), both on board the ISO satelHte. Although 
getting theoretical identifications and r edshifts of such objec ts 
is not an easy task (Lilly et al. ( |1999| ), Barger et al. ( |1999| )), 
the sources seem to be the high-redshift analogs of the local 
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"luminous" and "ultraluminous infrared galaxies" (LIRGs and 
ULIRGs). However, there are uncertainties on the properties of 
high-redshift dust, and on the mechanism that is responsible for 
dust heating (starbursts versus active galactic nuclei). In local 
LIRGs and ULIRGs samples, starbursts dominate, except in the 



brightest sources, with i > 3 x lO^^L© (Genzel et al. (1998) 



Lutz et al. (1998)). We will hereafter assume that starbursts are 



the component that dominates dust heating and we will neglect 
the influence of the AGNs. 

Optical studies are also showing that high-z galaxies se- 
lected from their optical prop erties, and in particula r the L BGs 
at z = 3 and 4 (Steidel et al. ( |l996t ), Madau et al. ( |1996| )), are 
so heavily extinguished that probing their properties requires 
multi-wavelength observations. Current estimates (Flores et al. 
( |1999[ ), Steidel et al. ( |1999[ ), Meurer et al. ( |1999[ )) give extinc- 
tion factors 0.1 < E{B — V) < 0.4, with a large scatter and the 
trend that the brightest objects are also the more extinguished 
ones. The net effect is that the SFR estimated from UV fluxes 
are too low by a factor 3 - 5 on average. Observations seem to 
show that dust is present at redshifts beyond 4 since it is seen, 
for instance, by its extinction of a gravitationally lensed galaxy 
at z=4.92 (Soifer et al. (1998)) and by its mm emission in radio 



galaxies and quasars (e.g. BR1202-0725 at z=4.69, McMahon 



etal. (1994)) 



So there is clearly a need for models of SEDs that take 
the extinction and emission effects of dust into account in a 



consistent way. Since the pioneering work by Tinsley (1972), 
Searle et al. (1973) and Huchra (1977), where the photometric 
evolution was restricted to the visible wavelength range, stel- 
lar population synthesis models were designed (Bruzual, 1983) 
and extended to include the nebular component and extinc- 
tion (3uiderdoni and Rocca-Volmerange, 1987). These mod- 



els were driven by the need to reproduce observed spectropho- 
tometric properties of local as well as more distant galaxies 
from the far-UV to the near-IR wavelength range. Spectropho- 
tometric evolution models were then improved (Bruzual and 
Chariot, 1993; Fioc and Rocca-Volmerange, 1997) to include 
tracks with various metallicities, and late stages of stellar evo- 
lution. They also g ained a better time resoluti on thanks to the 
isochrone scheme (Chariot and Bruzual, 1991 ). Work was also 



achieved to coup le chemical and spectrophotometric models 
self-consistently ( Arimoto, 1989 ). 

Other models were designed, where, in addition to chem- 
ical evolution, the photometric evolution of stellar popula- 
tions which heat the dust is taken into account (Franceschini 
et al. ^1991| ), ( |l994i ) and Fall, Chariot and Pei ( |1996| ), Silva 
et al. ( |1998 )). They all assume a simple relationship between 
the dust content and the heavy-element abundance of the gas 
(for details on the complex processes of dust grain creation 



and destruction, see e.g. Dwek (1998) and references therein) 



and do not have a detailed modelling of the re-processing of 
starlight by dust grains. On the other hand, models restricted 
to galaxy emission at IR/submm wavelengths have been devel- 
oped (see eg. Rowan-Robinson and Crawford (|1989|), Maffei 



( 1994 ), Guiderdoni et al. ( 1998| )) but none of them describe in a 
self-consistent and detailed way the stellar source that powers 



the emission. We would like to emphasize the need for mod- 
els which explicitly connect both UV/optical and IR/submm 
windows, because they are the only ones capable of correctly 
handling the luminosity budget of galaxies and hence their star 
formation histories. As a matter of fact, these two windows are 
naturally linked together For it is dust, which partially absorbs 
starlight coming out in the UV/optical window and re-emits it 
in the far-IR. 

In this paper, we present such a model, STARDUST, 
which contains up-to-date spectrophotometric modelling self- 
consistently coupled to chemical evolution and dust absorp- 
tion (following the same guideline as Guiderdoni and Rocca- 



Volmerange (1987), Franceschini et al. (1994) and Guiderdoni 



et al. (1998)). We further describe a new dust emission model 



following the lines of a model developed by Maffei ( 1 994 ) . Our 



approach differs from the one of Silva et al. (1998) in the sense 



that we empirically adjust our SEDs in the IR in order to match 
average IRAS colors, instead of trying to derive them from ba- 
sic physical arguments. What is lost in the understanding of 
the complicated physical processes involving dust is gained in 
the smaller number of parameters required by the model. Bear- 
ing this in mind, the model allows one to make consistent pre- 
dictions for the SEDs of galaxies in a very broad wavelength 
range. 

More specifically, most of the free parameters necessary 
to obtain such SEDs can be self-consistently computed in the 
explicit cosmological framework of semi-analytic models of 
galaxy formation and evolution (SAMs). However, we d efer to 
a companion paper ( pevriendt and Guiderdoni, 1999 ) a dis- 
cussion of the results obtained by implementing our SEDs 
in such SAMs. In Section 2, we briefly describe the spectro- 
photometric model. Section 3 deals with chemical evolution. 
Section 4 with the dust. In Section 5 we put the pieces together 
to build the complete spectral energy distribution of the galax- 
ies. Section 6 sketches out what a "dusty" high-redshift uni- 
verse might look like, and we eventually draw conclusions in 
Section 7. 



2. Computing Spectra of the Stellar Component 

2.1. Numerical Scheme 

The goal is to compute the stellar F^{t) flux contribution to 
the galactic spectrum at time t. Such a flux can be written as 
follows: 



mt) 



V'* {t — T)(p{m)f\ (m, T)dm dr 



(1) 



where rrid is the minimal star mass, to„ is the upper mass cut- 
off (respectively O.IM0, and 120 Af0 here), ip^,{t — r) is the 
number of Mq of gas that gets turned into stars at time (i ~ r) 
per time unit, (/)(m) ex m^^ (in the following we take x = 1.35 
for m in [nid, my] ( [Salpeter, 1955 )) is the IMF normalized to 
1 Mq, and f\{m, r) is the flux of a star of initial mass m at 
wavelength A and age r (t = corresponds to the zero age 
main sequence, and f\{in, r) = if t > t^,{ra) where t^,{m) 
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is the life time of a star of mass m). We shall neglect the nebular 
component in the following. 

The straightforward solution consists in discretizing this 
double integral, but leads to oscillations in the resulting flux 
mainly due to the rapidly evolving stellar stages (Chariot and 
Bruzual, 1991). The natural method to get around this, the 
so-called "isomass scheme" only discretizes the integral over 



mass ( GuiderdoniandRocca-Volmerange, 1987 ). However, the 
problem with such an approach is that it requires a very fine 
mass grid in order to achieve a good overlap of equivalent evo- 
lutionary stages for con secutive ma sses. Therefore, it is com- 
puter time consuming ( |^ioc, 1997 ). The other solution that 



leads to the same results but is much more computationally 
tractable is called the "isochrone scheme" and is based upon the 
discretization of the integral over time (Chariot and Bruzual, 
1991). It has the advantage of being directly comparable to 
color-magnitude diagrams of star clusters. The isochrones are 
of course carefully picked in order not to miss any stage of 
stellar evolution. We use such a scheme in the model presented 
here. 



Hence more recent grids of models, based on Geneva tracks 
and covering the evolution of low mass stars (0.8 to 1.7 Mq) 
from the zero age main sequence up to the end of the early- 
AGB, are included for Z = 0.001 and Z = 0.02 (Charbon- 
nel et al., 1996). For the late stages of other metallicities (hor- 
izontal branch, earlt AGB), we either interpolate or extrapo- 
late log Lboi, log Toff and log t versus log Z from the available 
tracks (Meynet private communication). The final stages of the 
stellar evolution (thermally pulsing-AGB and post-AGB) are 
not included. 

The main motivation for choosing the Geneva Group tracks 
library is the possibility to use Z-dependent stellar yields 
(Maeder (1992), Maeder (1993)) which insures consistency 
when spectrophotometric and chemical evolutions are coupled. 
The initial helium content of the gas from which stars form 
is computed starting with an initial helium fraction Y = 0.24 
after primordial nucleosynthesis and assuming an enrichment 
rate compared to metals of ^ = 3 (so that F = 0.3 for solar 
metallicity). 
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Fig. 1. Spectral evolution of an instantaneous burst of star for- 
mation with solar metallicity. The IMF is Salpeter, and the flux 
level corresponds to a 1 Mq galaxy. 



2.2. Library of Stellar Tracks 

Our model is mainly based on the so-called "Geneva tracks" 



(Schaller et al. (1992), Schaerer et al. (1993a) , Charbonnel et 



al. ( |1993| ), Schaerer et al. (1993b|)) (Z 
0.02, 0.04 and 0.8 < 



0.001, 0.004, 0.008, 
< 120 ). For stars less mas- 



2.3. Stellar Spectra 

An obvious advantage of using theoretical instead of empirical 
stellar spectra, is that the 3-parameter space of the theoretical 
spectra (log Toff, logiboi (or logi?) and logg), is similar to 
the one of the HR diagram. This avoids the use of the transfor- 
mations of bolometric luminosity log Lboi to visual magnitude 
My, and of logTcff to spectral types and luminosity classes, 
which are rather uncertain for the hottest and coldest stars. 

Here we use the theoretical fluxes grid of Kurucz 
( 1992| ) (hereafter K92) which covers all metallicities from 
logZ/ZQ — +1.0 to logZ/Z© = — 5.0 and 61 temperatures 
from 3500 K to 50 000 K. Each spectrum spans a wavelength 
range between 90 A and 160 /im, with a mean resolution of 20 
A. For the coldest stars (K and M-type stars) with T < 3750 K, 
K92 models fail to reproduce the observed spectra. Therefore, 
we prefer to use models comin g from d i fferent so urces. F or M- 
giants, the Bessel et al. ( |l989| ),( |l991b| ),( |l991a| ), ( |l991c| ) (here- 
after BBSW) models have been used. The wavelength range 
of BBSW fluxes is restricted to 0.49-5 /im. We have added a 
blackbody tail to the red end part of the spectra. On the other 
hand, BBSW M-giants present a spurious spike around 5000 
A, which we correct for, using the following procedure: we re- 
place the BBSW fluxes downward of 0.6 /^m by a handful of 



Gunn and Stryker (1983) M giants fluxes. To do this, we re- 



M/Mq 
sive than 1.7AIq, these tracks stop at the Giant Branch tip 



compute the Gunn and Stryker spectra of known spectral type 
and V-K colour for the effective temperature of BBSW spec- 
tra. Then we add a blackbody tail redward of 0.7 /im (Worthey, 
1994). As noted by this author, Gunn and Stryker fluxes are not 
suitable downward from 0.36 /im. We thus attach a blackbody 
for shorter wavelengths. In a similar manner, the M dwarf se- 
quence of BBSW is adopted (with log g = 4.5-5) but instead 
of adding Gunn and Stryker spectra at the posi tion of t he spike , 
we use M dwarf models computed by Brett (1995a) (1995b). 
Blackbody tails are added to complete the flux library. 
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For both BBSW and Brett models, we rebin the spectra us- 
ing the same wavelength grid as in K92. Furthermore, BBSW 
and K92 grids cover overlapping metallicity ranges which we 
interpolate to the five metallicities imposed by the Geneva 
tracks. 

2.4. Post-Starburst Evolution 

With the material that has been described in the previous sec- 
tions, one can build the spectrum of what is called a simple stel- 
lar population (SSP). This simply means that at time t=0, one 
turns all the gas {\Mq here) into stars which are distributed 
over the mass range [md,rau] according to the Salpeter IMF 
and one lets them evolve passively. Such a time evolution is 
shown in Fig. p]. 




Comparing the two figures also allows one to check the well- 
known trend that the higher the metallicity of the initial gas, the 
fainter the B band magnitude for t>l Gyr 

If order to compare such spectra directly with color data, 
one just has to convolve the SED with filters used to observe 
real galaxies. Color and magnitude evolutions are obtained for 
instantaneous bursts of star formation. This is shown in Fig. |[ 
A typical application for such a model is the study of globular 
clusters. 



> 

I 
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Fig. 3. Magnitude and color evolution of an instantaneous burst 
with a Salpeter IMF. The soUd curve corresponds to solar, 
the dotted line to twice solar, short-dashes to 2/5 solar, long- 
dashes to 1/5 solar and the dot-dashed line to 1/20 solar metal- 
hcities. Magnitudes are normalized to 1 M©. 



Fig. 2. Same as Fig. nlbut for a ^Zq burst. 



What is noticeable at first sight is the huge decrease of the 
flux with time which is all the more important when the wave- 
length is short, and simply reflects the fact that the most lumi- 
nous stars are the most massive ones and also the first ones to 
die. The second point is that one can clearly see the red super 
giant stars (RSG) population building up in the infrared long- 
ward of 1 /xm between 0.01 and 0.06 Gyr. 

These RSG stars are very sensitive to the metallicity of the 
initial gas from which they were born, in the sense that the 
higher the metallicity of this gas, the higher the flux emitted 
in the near-IR by the RSGs. This is clearly illustrated in Fig. 
for which the stars formed from a metal- poor gas. At i = 0.01 
Gyr, the flux at wavelengths greater than 1 /im is about an order 
of magnitude lower than the flux for an identical population of 
stars that formed from initial gas with solar metallicity (Fig. |l|). 



2.5. Comparing Codes 

The next interesting step is to see how our code compares to 
other population synthesis codes. Readers interested in more 
detailed comparisons (although on older versions of similar 
codes) are referred to Chariot et al. ( 1996 ) and the reviews by 
Arimoto ( |199^ ) and Chariot ( |1996| ). In Fig.| and |, we show 
a comparison of the photometry predicted by three different 
codes for three different metallicities. For the three codes, a 
similar standard Salpeter IMF has been used to distribute the 
stars in an instantaneous burst of star formation that takes place 
at i=0. 

The PEGASE code ( ^ioc, 1997 ) uses the Padova tracks 
(Bressan et al., 1993) for stellar evolution, with an initial he- 



lium fraction Y 

the end of the early AGB phase. It also includes the thermally 



0.23 and an enrichment of -^ = 2.5, up to 
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pulsing AGB ( Groenewegen and De Jong, 1993), and the post- 
AGB stages from Schonberner (|1983|) and Blocker (|1995|). The 



spectral libraries are fromLejeune et al. (1997). The GISSEL 



code (version 1998) uses tracks from Padova up to the early 
AGB phase which are completed till the end of the white dwarf 
cooling sequence and includes the thermally-pulsing, the post- 
AGB stages, as well as carbon stars (Chariot, private commu- 
nication). The spectral libraries are the same as for PEGASE. 
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Fig. 4. Comparison of GISSEL (version 1998) (dotted line) 
(Bruzual and Chariot, 1993), PEGASE (dashed line) (Fioc and 
Rocca-Volmerange, 1997), and our STARDUST (solid line) 
for an instantaneous burst with solar metallicity (lowest curve), 
2/5 solar (middle curve), and 1/5 solar (upper curve) in the 
different wavelength bands indicated on the panels. For a better 
visualization the 1/5 solar and 2/5 solar metallicities have been 
shifted by 2 and 4 magnitudes respectively. 

As can be seen in the different panels of Fig. ^ and ||, mod- 
els give very similar results in all bands, from the far-UV to the 
near-IR, although they use different stellar tracks and spectra. 
Most of the time, the differences are minor (around 0. 1 magni- 
tude), especially for colors and metallicities close to the solar 
value. However there are some discrepancies for early stages 
of evolution (t < 3Myr), in the far-UV for t > IGyr, and in 
the near-IR (mainly the K band). 

As far as the early evolution is concerned, although we 
are in fair agreement (better than 0.3 magnitudes) with the 
PEGASE and GISSEL models for all metallicities and wave- 
bands, the uncertainties in the early evolutionary stages of mas- 
sive stars (essentially mass loss) are responsible for the larger 
flux difference between the models. 




Fig. 5. Same as Fig. |4J. 

In the far-UV, after about 1 Gyr, differences are due to the 
uncertainties in the contribution of intermediate to low mass 
stars to the flux at short wavelengths. This is not relevant as 
long as one does not try to analyze E/SO galaxy spectra in 
this wavelength range, because, as mentioned e.g. by Rocca- 
Volmerange and Guiderdoni (1987), the level of such a flux 



is negligible compared to the far-UV flux produced by even a 
small amount of star formation. 

In the near-IR, the discrepancies (more pronounced for low 
metallicity bursts where they can attain 0.5 magnitude ) are 
due to the modelling of the final stages of stellar evolution 
(thermally- pulsing AGB and post-AGB) that are taken into 
account in PEGASE and GISSEL but not in our model. The 
thermally- pulsing AGB stars are responsible for the excess of 
flux of PEGASE and GISSEL in the near-IR (K band) and the 
redder V-K color from about 100 Myr to 1 Gyr. 

We want to emphasize the point that the comparison we 
have done here is by no means extensive, and is likely to de- 
pend on the IMF that one picks. Nevertheless, the remarkably 
good agreement observed in this specific case between the three 
totally independent codes is fairly encouraging. 

3. Chemical Evolution 

As previously mentioned, we use the Z-dependent stellar yields 



of Maeder (1992) with moderate mass loss for high mass stars. 



For intermediate mass stars, we use yields given by Renzini and 



Voli (1981). We also include type la supernova whose mod- 



elling and yields respectively follow Ferrini et al. (1992) and 



Nomotoetal. (1995) 



Our chemical evolution model is a simple closed box model 
(see e.g. Tinsley (1972)), where we assume that the sum of the 
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mass of stars Af * (t) and gas Mg (t) remains constant through 
time: 



Mg{t) + M^t) = Mtot , Mtot = C^ 



(2) 



We do not make the assumption of instantaneous recycling 
but we track metals ejected by stars at each time step and let 
new stars form out of the enriched gas, assuming instantaneous 
mixing of the metals. In other words, we solve at each time step 
the following equation : 



d{ZgMg)it) 

dt 



-Zg{t)^,{t)+£z{t) 



(3) 



where Zg is the mass percentage of element Z in the gas, V'* [t) 
is the amount of gas consumed by star formation, i.e. the SFR 
and is taken to be: 



Mt) 



Mg{t) 






(4) 



where i* is referred to as the characteristic time scale for star 
formation. 

Finally, £z (t) is the ejection rate of element Z at time t: 



£z{t) = 



d(logTO)^*(t — t{m))(f){m) 



X {(m — w{m))Z{t — t(m)) + mYzim)} , 



(5) 



where rrit is the minimum possible star mass ( i.e. the mass of 
a star having a lifetime of t). 

The first term between curly braces, (m — w{m))Z{t — 
t„i), is the ejected mass of initially present metals ( i.e. the 
difference between the initial mass of the star and the remnant 
mass, times the initial metallicity), the {t ~ t{m)) argument 
being due to the fact that stars of mass m which die at time 
t were born at {t — t{m)). The second term in these braces, 
mYz{m), is the percentage of the initial star mass transformed 
into new elements Z (i.e. the stellar nucleosynthesis). Yz{m) is 
referred to as the stellar yield. 

The results of such a chemical evolution model for star for- 
mation histories with various i, are shown in the upper left 
panels of Fig. ^ where we plot the time evolution of the gas 
fraction g{t) defined as g{t) = Mg{t)/Mtot and the gas metal- 
Ucity Z{t). 

4. Computing Dust Spectra 

4.1. Absorption 

To estimate the stellar flux absorbed by the ISM in a galaxy, 
one first needs to compu te its optical depth. As in Gui derdon i 
and Ro cca-Volmerange ( 1987 ), and Franceschini et al. ( 1991 ), 
(1994), we will assume that the mean face-on (z-axis) optical 
depth of the average gaseous disk at wavelength A and time t 
is: 



-(^^H't 



Zgit) 



{Nnit)) 



2.1 1021 at cm- 



(6) 




0.01 0.1 1 

t (Gyr) 



Fig. 6. The left-hand side panel shows time evolution of the 
gas fraction, metallicity and optical depth for the closed box 
model respectively described by equations g, |] and |[ Cod- 
ing for the different curves is given in the top panel, t, stands 
for the characteristic time scale for star formation as defined 
in equation ^ The right-hand side panel represents snapshots 
of mean face-on (z-axis) optical depth as a function of wave- 
length for a galaxy with a star formation time scale of 0. 1 Gyr 
according to formula H. The different lines represent the differ- 
ent snapshots at the time mentioned in the panel. All the curves 
of both panels are computed with /h = 1- 



where the mean H column density (accounting for the presence 
of helium) is written: 

{Nuit)) ~ 6.8 10^^g{t)fu atoms cm~2 . 



As noted by Guiderdoni ( 1987 ), a galaxy with g ~ 20% 
has (A^h) — 1-4 10^1 atoms cm^^ ^vhich is in fair agreement 
with the observational value for late-type disks, so that /h — 1 
for normal spirals (corresponding to rg/r25 — 1.6, see Guider- 
doni and Rocca-Volmerange (1987) for details). This is the 
value we will adopt in the "standard" case, but we will allow 
ourselves the possibility to increase the H column density to 
model ULIRGs for instance. In order to do this, we will just 
multiply the column density by a "concentration" factor /h 
which is inversely proportional to the surface occupied by the 
gas and accounts for the fact that star formation is more concen- 
trated in starbursts, than in normal spiral disks. Of course this 
is not completely satisfactory, and (A'h) should be physically 
linked to the size and mass of the gaseous disk in order to get 
rid of /h. The way we estimate sizes and masses of gas disks 
will be detailed in a companion paper (Devriendt and Guider- 
doni, 1999). 
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In equation g, the extinction curve depends on the gas 
metallicity Zg (t) according to power-law interpolations based 
on the Solar Neighbourhood and the Large and Small Magel- 
lanic Clouds, with s = 1.35 for A < 2000 A and s = 1.6 
for A > 2000 A (see Guiderdoni and Rocca-Volmerange 
( |1987 ) for details). The extinction curve for solar metallicity 
(AaMv)zo is taken from Mathis et al. ( |1983[ )P|. 

Examples of the behavior of the average face-on optical 
depth as computed using equation g can be read from Fig. ^. 
We note that with /h =1, the disks are optically thin during 
most of their evolution. 

4.2. Geometry 

Once one has computed the average face-on optical thickness 
of a galaxy, one needs to assume a given geometric distribution 
for the ISM, in order to compute dust obscuration from dust 
extinction properties. 



As in Dwek and Varosi ( 1996 ), we model a galaxy as an 
oblate ellipsoid where absorbers (dust) and sources (stars) are 
homogeneously mixed. Assuming a homogeneous absorption 
of the medium results in an overestimate of the absorption of 
the U V pho tons and is discussed in more details in Devriendt 
et al. ( 1998 ). Following the former authors we assume that the 
"thickness-to-diameter" ratio of our disk-shaped galaxies is ~ 
0.02. 



For spherical galaxies, Lucy et al. ( 198S ) generalized the 
analytic formu la giving obscura tion as a function of optical 
depth T^^ (t) ( [Osterbrock, 1989 ) to the case where scattering 
is taken into account via the dust albedo io\. Using Monte- 
Carlo simulations, Dwek and Varosi were able to show that 
the same formula is still valid provided one takes t\[€) ~ 
0.193 T^P^(i) ~ 2.619 r|(t) as the effective optical thickness 
of the disks at wavelength A and time t. The internal dust ob- 
scuration, i.e. the effective extinction, (averaged over inclina- 
tion angle i) is then given by: 



a\ 



1 - CJA + '^\a\ 



(8) 



(^aW), = -2.51ogio 

where 

' For wavelengths shorter than 912 A a second term should be 
added to account for h ydrogen absorption. This term is of the type 



(Devriendt etal., 1998) 



(iVH(t)>aH 



912. 



6(912. 



(7) 



is the Heavyside function, and cth = 6.3 x 10~^* cm'^ is the hydro- 
gen ionization cross section at the threshold. This term represents the 
internal extinction of the H ionizing continuum due to the presence of 
neutral hydrogen in galaxies which has to be added to the external ex- 
tinction of line-of-sight photons due to intervening gas rich systems. 
As it is dependent of the dumpiness of the gas and represen ts a small 
fraction of the reprocessed photons ( pevriendt et al., 1998 ), we will 
not consider this extra absorption term in what follows. 
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Fig. 7. Time evolution of the obscuration for a galaxy with dif- 
ferent geometries for the distribution of dust and stars. The star 
formation time scale is 0. 1 Gyr For panels on the left hand side 
/h — 1, whereas on the right hand side, /h = 100 (see text for 
a definition of /h). The thick dot-dashed curve is the polyno- 
mial fit given by Calzetti et al. (1994), rescaled to our Q(5500). 

Figure ^ shows the comparison of lyA\{t))i in our model 
with the commonly-used "slab" geometry. In a slab geometry, 
stars and dust are just distributed in the same infinite plane layer 
with the same vertical scale, so that lyA\(t))i is just given by: 



(^aW). 



-2.51ogi 



1 — exp(— -yi — u:\t\I COS %) 
vT^^73at^ / cos i 



where the albedo is phenomenologically taken into account 



through \j\ — LUX (Guiderdoni and Rocca-Volmerange, 1987). 
We also computed the absorption in a "screen" geometry, 
where the dust layer is in front of the stars, and in the "sand- 
wich" one, where the dust layer is sandwiched between two 
star layers. We do not show results for these two alternative 



cases, because, as shown by Franceschini & Andreani (1995) 



and Andreani & Franceschini (1996 ) for a sample of LIGs and 



normal spirals, these geometries result in too much and not 
enough starlight absorption respectively. On the other hand, the 
slab geometry (and the oblate ellipsoidal one which we use as 
our standard in what follows) seems to yield absorptions which 
are more in accordance with the data. 
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As can be noticed in figure ^ the results derived with the 
oblate ellipsoid model are quite similar to the ones obtained 
from a slab model. For a concentration factor /h = 1, i-e. for 
the so called "standard" model which corresponds to the av- 
erage obscuration of normal spiral disks, the agreement to the 



observational fit of Calzetti et al. (1994) is fairly good in the 



optical/near-IR wavelength range (4000 to 8000 A). But both 
models overestimate the UV extinction by quite a large factor ( 
about 1 magnitude at 1500 A), as can be seen on the left hand 
side of Fig. 0. This is due to the fact that the sh ape of the ob- 
scuration curve observed by Calzetti et al. ( 1994 ) is not a direct 
measurement of the extinction in a galaxy, but includes geomet- 
rical effects. It just represents the relative obscurations (obscu- 
rations at other wavelengths are compared to the obscuration 
at 5500 A). For different geometries, these relative obscura- 
tions scale differently with the optical depth. For instance, in 
the school case of the sandwich model only a maximum value 
of half of the emitted light can be absorbed, corresponding to 
a maximum value {A\{t))i = 0.75. Therefore, if absorption is 
already maximal at 5500 A, half of the light emitted at 2000 A 
will also be absorbed, even though the optical depth is much 
larger at this wavelength! To sum things up, the global effect 
is very much geometry dependent, and if one picks the "right" 
geometry for the distribution of dust and stars, one can either 
reduce the relative obscurations and smooth out any feature on 
the absorption curve ("sandwich like" effect), or increase the 
relative obscurations and enhance the features ("screen like" 
effect). 

This effect is illustrated on the right-hand panels of Fig. 0, 
where the same models of absorption are run, but with /h = 
100 which simply means that the optical depth is multiplied by 
100 at all wavelengths (see eq. g). One can see in Fig. ^ that 
the "slab" geometry exhibits a clear proclivity to reduce the 
relative extinction between 5500 A and shorter wavelengths, 
thus leading to a natural smearing of the extinction curve fea- 
tures (carbon bump at 2000 A) and to a better agreement (on 
the entire wavelength range) with the fit given by Calzetti et al. 
( 1994) for a sample of obscured starbursts. On the other hand. 



for the oblate ellipsoid geometry, the relative absorption in the 
same range remains quite insensitive to the value of the ab- 
solute optical depth. Considering our crude knowledge of the 
extinction curve of even close-by optically thin objects, and 
the sensitivity of absorption to geometry illustrated above, we 
adopt the oblate ellipsoid geometry in what follows, keeping 
in mind that it might overestimate absorption at wavelengths 
shorter than 3000 A. 

4.3. Dust Model 

In the previous section we have shown how to compute the 
total stellar luminosity absorbed by the ISM. Next step is to 
re-distribute this luminosity in the IR/submm window. 

In order to be able to compute such an emission spectrum, 
what one needs is to model the dust in the ISM. Our multi- 
component d ust mo del is mainly based upon the model by 
Desert et al. (1990). This model includes contributions from 



polycyclic aromatic hydrocarbons (PAHs), very small grains 
(VSGs) and big grains (BGs). 

Because of their small size (< 1 nm), PAH molecules are 
out of thermal equilibrium when excited by a UV/visible radia- 
tion field. Their temperature fluctuates and can reach a value 
well above the equilibrium temperature, naturally producing 
the bands at 3.3, 6.2, 7.7, 8.6 and 11.3 /im. Both VSGs and 
BGs are made of carbon and silicates. The main difference be- 
tween these grains is a size difference, though VSGs are prob- 
ably dominated by carbon whereas BGs seem to be dominated 
by silicates. VSGs have sizes between 1 and 10 nm. As a con- 
sequence, along with PAH molecules, they never reach ther- 
mal equilibrium. Therefore, their emission spectrum is much 
broader than a modified black body spectrum at a single equi- 
librium temperature. On the other hand, BGs, which have sizes 
between 10 nm and 0.1 /im (almost) reach thermal equilib- 
rium and can be reasonably described by a modified black body 
€^B^{T-qq) with emissivity e^, ex t/™ ( where the index m is 
such that 1 < 771 < 2) and temperature Tbg- Hereafter we take 

777 = 2. 



In contrast with Maffei (1994), we assume that the popula- 
tion of BG is divided in two components: 

- a cold component, with a fixed temperature Tbg = 17 K 

- a "starburst" component which has the same shape as the 
modified black body described above, but a higher Tbg to 
account for the fact that BGs receive extra heating from the 
radiation field in the star forming region. This component 
finds its justification in ob servations of e.g. the typical lo- 
cal starburst galaxy M82 ( Hughes et al., 1994 ), where the 
shape of the far-IR spectrum can be explained by a higher 
temperature (up to 40 K instead of 17 K) of BGs in thermal 
equilibrium with the radiation field. 



Contrary to Maffei (1994), we do not take into account the 
possible fluorescence of PAH molecules at wavelengths shorter 
than 3 microns. 



4.4. Assembling the Dust Puzzle 

The submm emission of galaxies, very sensitive to the spec- 
tral characteristics of dust, is more poorly known than the FIR. 
Submm fluxes of only a few tens of galaxies are available 
throughout the literature, and estimates of the amount of en- 
ergy released in this range are strongly discrepant (e.g. Chini 
et al. (|1986|); Stark et al. (|1989|); Bales et al. (|1989|); Chini and 



Ki-ugel (1993)). However, Stickel et al. (1998) have shown that 



data from the ISOPHOT shallow survey at 175 /itm are incon- 
sistent with a large population of cold galaxies undetected by 
IRAS. 

Guided by all these data, one can compute emission spec- 
tra of galaxies over the whole wavelength range (near-IR to ra- 
dio). We do this by adding the different contributions from the 
populations of dust grains (see previous subsection) to the cor- 
responding post-absorption synthetic stellar spectra. The dust 
dominated part of the spectra is picked among the spectral li- 
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brary (see previous subsection) following a method similar to 
the one developed by Maffei (|1994|). 




Fig. 8. Spectral energy distribution in the IR window for galax- 
ies with different total infrared luminosities ranging from 
IO^Lq (bottom spectrum) to lO^^L© (top spectrum). Solid 
line is our model, and the dotted curve is the model of Maf- 
fei (1994). 



and more mid-IR to submm emission due to the VSGs and 
BGs. 

To further extend the wavelength range, synchrotron ra- 
diation is added to the spectra. This non-thermal emission is 
strongly correlated with stellar activity and, as a consequence, 
with IR luminosity (see e.g. Helou et al. (1985 )). According to 
observations at 1.4 GHz, this correlation is L^{\A GHz) = 
Xir/(3.75 10^2 X 109), where im is in W, i^go = 3.75 lO^^ 
Hz is the frequency at 80 jim. and q ~ 2.16 is determined from 
observations. Then we assume that we can extrapolate from 21 
cm down to ~ 1 mm with a single average slope 0.7, so that 
L^ = L^{IA GHz){y/lA GHz)-^'^. 

This is the way we proceed for "normal" infrared galaxies 
as well as for the models that are presented here unless the con- 
trary is mentioned. For the sample of ULIRGs of Rigopoulou 



et al. (1996) we notice that the average slope is somewhat shal- 
lower, around 0.46 and the parameter q is more like q ~ 2.85. 
As the sample is not large enough for this correlation to be 
firmly established, we fit each ULIRG individually with spe- 
cific q parameter and slope, in order to reproduce the observed 
radio flux. We also note that a steeper slope for some objects in 
the sample (e.g. Mrk 273) probably indicates a non-negligible 
contribution in the radio flux coming from an active nucleus. 
Nevertheless, we assume that the dominant contribution in all 
the galaxies of this sample comes from the star forming region, 
thus neglecting a possible pollution by the active nucleus. This 
indeed seems to be a good approximation, in view of results 
derive d from a recent study of a sample of 60 ULIRGs by Lutz 
et al. (1998), as we expect 80 % of the objects to be powered 
by star formation in the IR. 



This method consists in using the observational correla- 
tions of the IRAS flux ratios 12/im/60/im, 25/im/60/im and 



lOO/im/60/^m with Ljr ( Soifer and Neugebauer, 1991 ) to de 



termine the contribution of each component of the dust model. 
These correlations are extended to low Lir w ith the samples of 
Smith et al. ( |l987| ) and especially Rice et al. ( |l988| ). We would 
like the reader to bear in mind that these samples are rather 
small, so that the colour-luminosity correlation is somewhat 
loos er At submm wavelengths, the samples of R igopoulou et 
al. (1996) and Andreani and Franceschini (1996) are used to 



derive a correlation of the 350/im/60/im flux ratio with Lir. 
As a result, we have now four color ratios which we use to 
caUbrate our four components. 

In Fig. [sj we show the SEDs predicted by our model as a 
function of total infrared luminosity. First of all, we notice an 
overall robustness of the derived SEDs for ijR > IQ^Lq, in 
the sense that the shapes and emission maxima of our SEDs 
are very similar to those first derived by Maffei (1994). We 



have slightly more flux at wavelengths between 200 and 2000 
/xm for galaxies with Lir < lO^^i© which is the effect of 
our different modelling of the BGs. As a matter of fact, the 
cold part of the BG component is responsible for the excess of 
flux (about a factor ~ 2 at 850 ^m) that is seen in Fig. |8) For 
galaxies with Ljr < lO^L© we also have less PAH emission 



5. Building tlie Full Synthetic Spectra 

Now that we have shown how to compute both the stellar spec- 
tra and the dust spectra and how to link these two components 
in a self-consistent way, we just piece them together in order to 
obtain the full galaxy spectra from the far-UV to radio wave- 
bands. This is shown in figures |^, [I^ and |ri|, for different types 
of galaxies (normal spirals and starbursts with normal and high 
extinction) and for different kinds of geometries. In these fig- 
ures, we also overplot the original spectra as they come out 
from the spectral synthesis models without extinction. 

The important thing to notice is that for normal spirals, the 
amount of luminosity absorbed by dust is about half the lu- 
minosity produced by star formation, whereas for a starburst 
where the optical depth is calculated with the same prescrip- 
tion, the luminosity released in the IR is already twice as large 
as the luminosity emitted in the optical. This is mainly due to 
the fact that, in starbursts, star formation takes place on a time 
scale which is much shorter than in normal spirals, so that the 
stellar population has less time to grow old and it emits most of 
its flux in the UV/optical window where absorption is higher. 
The other important point to notice is that, when the objects 
remain marginally optically thin (with t| < 1), which is the 
case for the objects plotted in Fig. ^ and |l^ the UV/near-IR 
SEDs (except maybe at wavelengths < 2000 A) are fairly ro- 
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Fig. 9. Snapshot of the full wavelength range synthetic spec- 
tram of a typical spiral galaxy with a star formation time scale 
of 3 Gyr and /h = 1 taken at time f = 10 Gyr, correspond- 
ing to Ty ~ 1 (see bottom left panel of Fig. ^). The ini- 
tial gas mass available for star formation has been renormal- 
ized to A/g(0) — Mtot — 1 ^Iq (see eq Q), the real initial 
value (which has a direct influence on the IR spectrum) being 
Afg(O) == 10^° Mq. The different curves represent different 
geometries of dust and stars. Thick solid line is the homoge- 
neous oblate ellipsoid mix, and thick dot-dashed line assumes 
a screen geometry. As a guideline, we also plot (thin dashed 
curve) the spectrum without any absorption. 



bust as they do not depend sensitively on the geometry of the 
distribution of stars and dust. 

For compact starbursts, the optical depth can increase by 
about two dex (we model these objects by putting our con- 
centration factor /h = 100), and, as shown in Fig. [ri|, two 
important things happen. First of all, the ratio of the far-lR 
to optical fluxes is multiplied by a hundred, which was to be 
expected. Second, and less trivial is that, depending on the ge- 
ometry adopted, the shape of the optical spectrum can undergo 
drastic changes. For instance, if one compares figures |l^ and 
Ml it is obvious that, for an oblate ellipsoid geometry, the shape 
of the stellar spectrum is quite independent of the value of the 
concentration parameter /h, and the optical SED just under- 
goes an attenuation of a couple of dex. On the contrary, for an 
extreme geometry like the screen model, the change is so dras- 
tic that the stellar spectrum that was present with /h ~ 1, has 
totally disappeared for /h — 100. Of course, this is purely a test 
case, for which the radiative transfer can be solved analytically, 
and the absorption (neglecting inclination) varies exponentially 
with the optical depth, but nevertheless, it is very illustrative of 



Fig. 10. Snapshot of the full wavelength range synthetic spec- 
trum of a typical starburst galaxy with a star formation time 
scale of 0.1 Gyr and /h = 1 taken at time t = 10 Gyr and 
corresponding to Ty ~ 1 (see bottom left panel of Fig. ^. The 
initial gas mass available for star formation has been renor- 
malized to Mg(0) = Mtot = 1 Mq (see eqg), the real initial 
value being Mg(0) = 10^^ M©. The line coding is the same 
as in Fig. p. 



the sensitivity of the optical spectrum to the dust distribution 
relative to stars. The explanation of this dependency is that in 
the case of the oblate ellipsoid geometry, we have a homoge- 
neous mix of dust and stars, so when the optical depth increases 
we still see the light coming from the stars contained in the 
"outer peel" of the galaxy. On the contrary, for the screen ge- 
ometry, the wall of dust in front of the stars just thickens, and 
does not let any light go through at all. 

Another point we would like to emphasize is the robust- 
ness of the dust emission spectrum. As a matter of fact, the 
level of flux for all wavelengths > 10/im differs by less than 
20% between "compact" starbursts (Fig. |l^ and "normal" star- 
bursts (Fig. |n]) even if the absorption of the starlight is two 
dex higher for a compact starburst. This is due to the fact that 
the starlight is absorbed in a narrow wavelength range (typi- 
cally 3 /im wide) and re-emitted over a much larger wavelength 
range (typically 1000 /im wide). Moreover, the extra quantity 
of starlight absorbed by the compact starburst represents a very 
small fraction of the bolometric luminosity Lboi of the object 
(and hence a very small fraction of the total Ljr ~ 2/3Lboi)- 
These combined effects have the dramatic consequence that 
one cannot predict the visible spectrum of an ULIRG! How- 
ever, one should be aware that the shape of the IR/submm SED 
is very likely to be altered in compact starbursts, due to absorp- 
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tion of dust re-processed photons by other dust grains. When 
the optical depth reaches values as high as w 100 in the V band, 
one can expect it to still be of the order of 10 at 3 /im, and 
more or less of the same value up to 30 /im because of the 
silicate absorption features around 10 and 25 fim (see Fig. ph. 
Our simple model does not account for such self-absorption 
explicitly, although one could argue that in principle it is more 
or less naturally accounted for, simply because of the way we 
build our spectral library from the phenomenological colors of 
IRAS galaxies. To sum things up, this effect should result in 
an important decrease of the strength of PAH features in highly 
obscured starbursts, which emit more than 90 % of their lumi- 
nosity in the far-IR, and for this reason, it should be present in 
IRAS observations of ULIRGs which we use to build our spec- 
tral library. We therefore believe that our conclusion about the 
overall robustness of the IR SEDs still holds when dust self- 
absorption is explicitly taken into account. 



Fig. 11. Same spectrum as in Fig. |lO| but for /h = 100 thus 
corresponding to Ty ~ 100. The real initial gas mass available 
for star formation is still Mg{0) = Mtot = 10^^ Mq (see eq 
g). For the screen geometry (dot-dashed line), the level of the 
flux below 2 /im is negligible. 
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Fig. 12. Time evolution of the luminosities and colors of a 
galaxy with different star formation time scales and /h = 1- 
Coding for the lines is the same as in Fig. H top left panel. 
The initial gas mass available for star formation is Afg(O) — 
Mtot = 10^^ Mq renormalized to 1 Mq, and Lq stands for 
solar bolometric luminosity. 
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Fig. 13. Same as in Fig. hjbut with /h = 100. 



In Fig. |I^ and |l3|, we quantify these changes a bit further 
by plotting integrated quantities and broad band colors. As ex- 
pected, before 3 Myr, there is no IR flux, because no stars have 
yet exploded. So the ISM is not enriched with metals yet, and 
there is no dust. After 3 Myr, the metallicity starts growing, 
along with the dust content and the IR luminosity of the galaxy, 
while in the same time more gas gets turned into stars. As the 
optical depth is a function of both the metallicity and gas con- 
tent, the IR luminosity reaches a maximum and then decreases. 
In fact, things are more complicated because this behavior de- 
pends on the stellar emission spectrum. But, as shown in the 
same plot (Fig. |l2|), the optical luminosity has a similar evolu- 
tion, which is the result of a competition between stars forming 
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at a rate which decreases with time because it is proportional 
to the gas content of the galaxy, and the ageing of these same 
stars. 

What one notices at first glance when one compares Fig. |l2| 
and |l3|, is that the global behavior is similar to what has previ- 
ously been mentioned. The shapes of the spectra are very much 
alike (mainly because the oblate ellipsoid geometry does not 
alter significantly the features and the extinction curve as pre- 
viously discussed) in the sense that the colors (especially B-V) 
are not noticeably different. Otherwise, the variations of the 
plotted variables are those expected when one boosts extinc- 
tion: the IR luminosity increases, the optical luminosity drops, 
and the absolute magnitudes in B and K increase. 

The next step after checking the self-consistency of the 
model consists in testing the spectra against real data in the 
full wavelength range. 
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Fig. 14. Flux ratio IRX= logig (FIR/F220W) as a function of 
UV spectral slope /?. Data is from Meurer et al. (1997) and 
Kinney et al. (1994). Galaxies with different star formation 
time scales are overplotted. Coding for the lines on top panel 
is : solid for i* = 0.1 Gyr (thin line /h = 1 and thick line 
/h = 100), dotted line for i* = 1.0 Gyr, and dashes for 
i* = 10.0 Gyr (the thickness coding remaining the same). For 
the bottom panel, isochrones are plotted for different t* rang- 
ing from 0.1 to 10 Gyr. Coding for times at which isochrones 
are taken is indicated in the panel. Once again a thin line stands 
for /h = 1 and a thick line for /h = 100. See text for more 
details and figure |I5| to estimate the starting values of f3 and 
IRX in the bottom panel. 

Fig. |l^ and |lj show how the models compare to observa- 
tions in the UV and far-lR simultaneously. The InfraRed eX- 



cess (IRX) — which is just the decimal logarithm of the ratio 
of the far-lR flux (as estimated from IRAS observations) and 
the flux measured in the F220W filter of the Hubble Space Tele- 
scope (see Meurer et al. (1995)) — , is shown both as a function 



of time (Fig. h5b and as a function of the UV spectral slope /3 
— which is defined as the power law index of the UV contin- 
uum Fa (A) oc A'' between 1250 A and 2600 A , (see Calzetti 
et al. ( |l994l ) for details) — (Fig. ^. 
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Fig. 15. Time evolution of the flux ratio IRX and UV spectral 
slope f3 of galaxies with different star formation time scales. 
Coding for the lines on both panels is as follows: solid line 
stands for t, — 0.1 Gyr, dotted line for i* = 1.0 Gyr and 
dashes are for t^ — 10.0 Gyr. In the top panel, the thin line 
represents models without extinction, the thick line, models 
with /h = 1 and the very thick line, models with /h ~ 100. 
Obviously, for the bottom panel, there is no thin curve, but the 
thickness coding remains the same in the other two cases. 



There are a few important things to notice on these figures. 
In particular, (upper panel of Fig. |l5|), one can see that, for 
our geometry (oblate ellipsoid), the UV spectral slope remains 
negative whatever the value of /h, and that, in order to get (al- 
most) flat spectra in the UV (i.e. /3 ~ 0) for starbursts, we need 
to invoke a higher extinction by setting /h = 100. 

It is striking that it is almost impossible to obtain positive 
values of f3 with reasonable values of the star formation rate 
and extinction. Of course, this statement depends on the ge- 
ometry as well as on the IMF, but one would have to come 
up with a very ad-hoc combination of these two parameters in 
order to change the sign of (3. Furthermore, as shown on the 
same figure (bottom panel), higher values of (3 are not corre- 
lated with a strong IRX. They are reached before the IRX at- 
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tains its maximum level. This just means that if one demands 
that the geometry be the homogeneous oblate ellipsoid and the 
IMF be Salpeter, one would have to advocate a column density 
at least a 1000 times larger in starbursts than in normal spirals 
(the metallicity dependence becomes second order due to the 
fact that the peak in the time evolution of f3 appears very early) 
in order for the extinction to be sufficient to change the sign of 
/3 (/h = 1000 in this case, and (3 reaches a maximum value of 
0.27). 

However, as can be read in Fig. |l^ the vast majority of 
galaxies have negative UV continuum slopes, and therefore, 
the models are able to span the whole data range quite natu- 
rally. As a matter of fact, the spectral slope of the extinction- 
free starburst is /? ~ —2.8 for t < 0.02 Gyr , in agreement with 



what is generically used (Steidel et al., 1999). Furthermore, the 
predicted reddening is 0.04 < E{B — V) < 0.4 (mainly de- 
pending on /h) for t <0.1 Gyr, in fair agreement with what is 
derived in LBGs under the assumption of a universal (averaged 
over 0.1 Gyr) /? ~ -2.5 (|Steidel et al.,"T999|). 



The results of the method we use to obtain the fit can be 
summarized in the following way, depending on the nature of 
the fitted objects: 

- In the case of the ULIRG sample (in which we include M82 
which is not an ULIRG), we find that for most objects, an 
old population (with a typical age of 10 Gyrs) with the 
same extinction as the starburst population (/h — 100) is 
required to account for the excess of flux in the near-IR 
(around 1-2 /im). This is due to the fact that even in our 
simple modelling, we cannot assume that the measured 
flux is entirely due to the burst, and we have to take into 
account the previous star formation history of the galaxy. 
Clearly, in some cases, there is a degeneracy between i* 
and this population. This can be understood if one keeps 
in mind that increasing i* simply means that the stellar 
population builds up more slowly, and consequently has 
time to grow old and therefore emits more flux in the 
near-IR. 



5.7. Fitting Observed Galaxies 

The next step is to test the ability of our models to describe 
multi-wavelength flux measurements of a broad class of ob- 
jects, from normal galaxies that emit approximatively 30 % of 
their bolometric luminosity in the IR, to monsters which have 
Lir/Xb ^ 100. 

To achieve this comparison, we simply select from the ht- 
erature a sample of galaxies for which fluxes in the UV, optical, 
IR and submm are available. We also pick the sample in order 
to have galaxies with total IR luminosities (as estimated from 
IRAS observations, see e.g. Sanders and Mirabel (1996)). All 
the galaxies we show here are taken from samples described in 



Rigopoulou et al. (1996) for ULIRGs, and Boselli et al. (1998) 



for normal spirals. We also add to the list the close-by LIRG 
M82 for which our multi-wavelength data availability criterion 



is fulfilled (see Hughes et al. (1994 ) and references therein). We 
sum up the properties of the selected objects in table |l} 



The parameters that we consider free when we perform 
the fit are the star formation characteristic time-scale t*, and 
the age of the galaxy t. The fits are obtained by a standard 
X^ minimization procedure, although we are aware that it has 
not much statistical significance, mainly because the number 
of data points used to perform the fit is too small (at most 17 
fluxes at different wavelengths). Moreover, the assumption that 
errors are distributed normally in the statistical sense (which 
is a requirement for a real x^ procedure to be valid) probably 
does not hold. Typical values of the parameters t and i» for 
the objects represented in the different panels of Fig. |l6|can be 
read in table [l|, along with typical properties derived from the 
best fit model. For the remaining objects of table |lj spectra are 
available from the authors upon request. 



- For the sample of normal spirals, we allow the same pa- 
rameters i* and t to vary, but we take /h = 1- We then add 
in a few cases a heavily extinguished starburst (/h = 100) 
to account for the missing far-IR flux. Once again, this is 
due to the fact that having just one characteristic time scale 
for star formation in our model implicitly assumes that, for 
these objects, all the measured flux comes from the quies- 
cent continuous star formation history of the galaxy, which 
might describe the average properties of a sample of galax- 
ies, but is too crude when one deals with individual objects 
in an IR-selected sample. 

From all the individual objects fitted, we were then able to 
extract a sequence of galaxies with different total IR luminosi- 
ties. This is shown in Fig. |l^ (in the spirit of Sanders & Mirabel 
(|19M)). 



One can see in this figure that, as the IR luminosity in- 
creases, the near-IR to UV spectral slope becomes shallower, 
and the characteristic features like the 4000 A break are 
smeared out. This is the result of a combination of lower ab- 
sorption and older stellar population for the sample of normal 
spirals as opposed to higher extinction and younger stellar pop- 
ulation for the ULIRGs. However, the figure is somewhat mis- 
leading, because one could also see a luminosity sequence in 
the UV/near-IR SEDs that correlates with Ljr . As pointed out 
earlier in the text, there is no such sequence, for the simple 
reason that the UV/near-IR SED is extremely sensitive to the 
relative distribution of stars and dust in heavily obscured ob- 
jects. For instance, Arp 220 has a SED in this wavelength range 
similar to M82 and its Lir is about 100 times larger (see ta- 
ble nl). Finally, the peak of the infrared emission gets shifted 
towards higher temperatures or equivalently towards shorter 
wavelengths as the total IR luminosity of the object increases. 
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Fig. 16. Best fit model (see text for details) for galaxies extracted from a sub-sample of the Virgo Cluster Catalogue (Boselli et 
al., 1998), and ULIRGs taken from the sample described in Rigopoulou et al. (1996). As mentioned in the text, this is illustrative 
of the ability of the model to capture the characteristic features of the objects and there is an important degeneracy between the 
different parameters (t* and t) of the fit. Objects are ordered with increasing Lir from top to bottom and from left to right. 



6. A Guide to High-z Galaxies 



make predictions for objects which look like these local galax- 
ies. 



More specifically, we show in Fig. 18 how the brightest 



The model described in the previous sections is designed to ULIRG from our sample would look at different redshifts. We 
matchavailablemulti-wavelengthdataof local galaxies. In this would like to draw the attention of the reader to the crucial 
section, we use the template spectra derived from the model to effect of the so-called negative k-correction at submm wave- 
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Table 1. Properties of the different objects for which individual fits were achieved. Values given for the absolute magnitudes 
are derived from the best fit spectrum, i.e. the spectrum of the galaxy which age t and characteristic star formation timescale 
i* minimize the ^^. For Ljr, values reported are obtained by integrating the spectra between 3 /im and 1000 /im, whereas, for 
Lqpt, we have integrated the spectra between 100 A and 3 /im. SFR and ry are indicative values that are derived from the best 
fit model. 



Object Name 


U (yr) 


i(yr) 


D (Mpc) 


Mb 


Mk 


Lm (Lq) 


LOPT (Lq) 


SFRCMeyr^") 


TV 


VCC 0836 


10'-* 


6.7 10^ 


17 


-14.7 


-18.9 


1.7 10** 


3.0 10** 


1.0 10-2 


0.6 


VCC 1003 


10« 


1.4 10^° 


17 


-14.2 


-18.2 


2.1 10' 


1.8 10** 


2.3 10"^ 


0.1 


VCC 1043 


10^ 


1.4 10^° 


17 


-14.5 


-18.5 


4.1 10' 


2.4 10** 


3.0 10"^ 


0.1 


VCC 1554 


lOi" 


3.2 10^ 


17 


-16.0 


-18.8 


6.9 10** 


5.2 10** 


1.0 10"^ 


0.6 


VCC 1690 


10^ 


9.3 10^ 


17 


-14.8 


-18.8 


7.8 10' 


3.1 10** 


5.8 10"^ 


0.3 


VCC 1727 


10« 


1.1 10^° 


17 


-14.5 


-18.6 


5.0 10' 


2.5 10** 


3.8 10-=* 


0.2 


VCC 1972 


10^" 


3.4 10^ 


17 


-15.4 


-18.4 


3.1 10** 


3.2 10** 


6.0 10-2 


0.6 


VCC 1987 


10^0 


3.3 10^ 


17 


-15.3 


-18.2 


2.7 10** 


2.9 10** 


5.4 10-2 


0.6 


M82 


5 10** 


9.0 lO'' 


4.2 


-19.3 


-22.9 


6.0 10^" 


1.6 10i° 


9.1 10° 


2.2 


Arp220 


10** 


5.0 10' 


108 


-19.4 


-23.7 


1.6 10^^ 


2.2 10^° 


2.7 10^ 


75 


Mrk 273 


10** 


3.0 10' 


222 


-20.7 


-25.0 


2.6 10^^ 


7.1 10^° 


4.2 10^ 


34 


UGC 05101 


10** 


2.0 10' 


240 


-20.9 


-25.3 


1.7 10" 


9.3 10^° 


3.0 102 


15 


IRAS 05189-2524 


10** 


3.0 10' 


256 


-21.3 


-26.1 


2.7 10" 


1.6 10" 


4.6 102 


34 


IRAS 08572+3915 


10** 


2.0 10' 


349 


-20.7 


-23.8 


2.4 10" 


4.7 10^° 


4.5 10^ 


15 


IRAS 12112+0305 


10** 


5.0 10' 


435 


-20.4 


-24.6 


3.4 10" 


5.1 10i° 


6.6 10^ 


56 


IRAS 14348-1447 


10* 


3.0 10' 


495 


-20.9 


-25.0 


4.0 10" 


7.8 10^° 


6.5 10^ 


34 


IRAS 15250+3609 


10* 


3.0 10' 


319 


-20.2 


-24.4 


2.0 10" 


4.3 10^° 


3.3 10^ 


34 



lengths which makes such objects as luminous at z ~ 5 as at 
z ~ 0.5 for wavelengths ranging between 400 and 900 fim. 

In Fig. p9| and |20|, we show optical/near-IR properties that 
various classes of local objects (from small normal spirals to 
bright ULIRGs) would have, if they were located at high red- 
shift. The two different pictures represent two different cos- 
mologies (Standard CDM and ACDM, see figure captions for 
details), and the differences induced by this are of the order of 
0.5 magnitudes in the various filters. Also the dot-dashed line 
shows the evolution effects one expects to measure on the spec- 
trum of a spiral galaxy which forms at redshift 5 with an initial 
gas mass Mg{0) = lO^o A/©. 

Fig. |T] and g^ ^r^ far-IR/submm properties of the same 
classes of objects. Here again, the important thing to notice is 
that the change in fluxes due to the different cosmologies is 
quite small (less than 30 %) at any wavelength. 

Eventually, Fig. ^ shows the time evolution of UV lumi- 

(left panel) and IR luminos- 



nosities X-jsnnL 



Aifinni> 



2800-t^A28oo' ^1600-t^Ai 

ity Lijx (right panel), considered as functions of the star forma- 
tion rate. We remind the reader that the SFR time evolution is 
driven by equation H, and that our galaxies are considered as 
isolated objects ("closed box" approximation). The interesting 
thing to note is the tightness of the correlation between iiR and 
SFR over the five orders of magnitudes spanned by the galaxies 
we picked. 

One could point out that there seems to be a correlation 
between the UV luminosities and SFR too, but, looking care- 
fully at Fig. ^ this holds only for SFR smaller than a few 
Mq yr-^. This is obviously due to the high extinction of all 
the objects that are massively forming stars (SFR > 100 7\/q 



yr~^): their UV luminosities are not significantly different from 
objects with more normal SFRs (up to 10 Mq yr"^). How- 
ever, we emphasize that this should be viewed as a theoreti- 
cal correlation, because the points plotted on Fig. ^ are de- 
duced from the best fit model. Nevertheless, the mathematical 
expression we get from simply fitting this "data" with a power 
law is LiR {Lq) = 7.7 x 10^ (SFR / Mq yr-i)0-95 
where the exponent of the SFR is remarkably close to 1. 
The same fitting method for the UV luminosities yields : 
A2800 iA.«oo (^o) = 2.0 X 108 (SFR / Mq yr-^-^^ and 
Aieoo iAieoo (io) = 2.4 X lO8(SFR/M0yr-i)O-58with 
exponents fairly different from unity. Furthermore the scatter 
is definitely more important. AU this points out that it is in- 
deed very important to have a fair estimate of the Ljr of high- 
redshift galaxies if one wants to derive realistic values for their 
SFRs. 

7. Conclusions 

In this paper, we have described a simple model which en- 
ables us to reproduce the far-UV to submm SEDs of observed 
galaxies, with a possible extension to radio wavelengths. The 
specificity of our model is the fact that it self-consistently links 
starlight and light reprocessed by dust to derive the global spec- 
tra of galaxies. We hereafter summarize the main assumptions 
of our work that attempts to be as phenomenological as possi- 
ble. 

To obtain SEDs of synthetic stellar populations, we couple 
spectrophotometric and chemical evolutions in a simple way, 
under the assumptions of a constant IMF (we take Salpeter 
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Fig. 17. Sequence of galaxies taken from our best fit models 
(see table |l]) with different total IR luminosities as in Sanders 
and Mirabel 1996. Distances to spiral galaxies in the Virgo 
Cluster are taken to be 17 Mpc. For the ULlRGs, distances are 
estimated by Sanders et al. (1988), rescaled to Ho = 50 km 
s~^ Mpc~^. Coding for the lines gives the name of each ob- 
ject: thin solid line: VCC 1003; thin short-dashed line: VCC 
1727; thin dotted hne: VCC 1554; thin long-dashed line: VCC 
1987; thick dot-dashed line: M82; thick short-dashed hne: 
Arp 220; thick dotted line: IRAS 15250+3609; and thick sohd 
line: IRAS 12112+0305. 



IMF) and no gas inflows/outflows (closed-box). Once stellar 
SEDs are computed, we use a phenomenological extinction 
curve that reproduces the observed metal-dependent trend in 
the SMC (Z = 1/5 Zq), LMC (Z = 1/3 Zq) and Milky Way 
(Z = Zq). The optical thickness can be computed e.g. un- 
der the assumption of no gas inflows/outflows. A simple oblate 
geometry where stars and dust are homogeneously mixed fi- 
nally gives the obscuration curve. The absorbed luminosity is 
redistributed at IR/submm wavelengths with a phenomenolog- 
ical model involving various dust components that reproduce 
the local correlation of IRAS and submm flux ratios with IR 
luminosity. The stellar and dust SEDs are then coupled. In this 
simple approach, the predictions depend on the SFR timescale, 
galaxy age, and size of gas disk that is simply parameterized 
with the fudge factor /h. 

We have shown that the results derived separately for the 
UV/optical part of the spectrum and the far-IR/submm lumi- 
nosity distribution are robust in the sense that they are very 
close to the results obtained independently by other authors. 
The predicted stellar SEDs are in good overall agreement with 
other predictions available in the literature, in spite of different 



Fig. 18. Observer- frame spectrum of best fit model of the 
ULIRG galaxy IRAS 14348-1447 at increasing redshifts, for 
a cosmology where Ho — 50 km s~^ Mpc~^, ilo — 1 and 



origins for the stellar data. The differences are mainly due to 
the fact that the late stages of stellar evolution (TP-AGB and 
post-AGB) are not included in our modellin g. The dust SEDs 
are in good agreement with those of Maffei ( 1994 ), who does 
not use submm data and takes slightly different components 
for the dust. The differences amount to at most a factor 2, and 
decrease with increasing luminosity. 

Furthermore, we have demonstrated that the results of the 
observational and theoretical w orks o n extinction by Calzetti 
et al. ( 1994 ) and Meurer et al. ( 1995 ) can be naturally inter- 
preted in the framework of our model as geometrical effects 
in the relative dust/star distributions. Indeed, in the optically- 
thick regime with optical depth Ty ~ 100, we observe a flat- 
tening and smearing of the {A\)i curve — as reported by 
Calzetti et al. (1994) — that depends on the way dust and 



stars are distributed. The obscuration curve is in good agree- 
ment with Calzetti's, though it is less "grey". A possible ex- 
planation is that Calzetti's starbursts are more metal-rich on 
average, than our model starbursts that start from zero metal- 
licity. The spectral slope of the extinction-free starburst is 
/3 ~ —2.8 for t < 0.02 Gyr , in agrement with what is gener- 



ically used (Steidel et al., 1999). The predicted reddening is 
0.04 < E{B -V) < 0.4 for t <0.1 Gyr, in fair agreement 
with what is derived in LBGs under the assumption of a uni- 
versal P ~ -2.5 ([Steidel et al., 1999|) for the UV slope of 



galactic spectra without extinction. The relation of the IRX to 
f3 reproduces the observed trend, provided /h = 100 for the 
most extinguished starburts, and points towards an age of i = 
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Fig. 19. Apparent magnitude for 3 objects, the "normal" spiral 
VCC 836 of the Virgo Cluster (sohd line), the close-by LIRG 
M82 (dotted hne), and the ULIRG IRAS 14348-1447 (dashed 
line), as a function of redshift, for an Einstein-de Sitter cos- 
mology where Ho = 50 km s^^ Mpc^^, Oq = 1 and Oa = 0. 
Symbols are galaxies from table |l|. Dot-dashed line is a model 
spiral with t^ — 3 Gyr and an initial gas mass available for star 
formation of Mg = 10^" M0 that formed at redshift 5, and 
for which the evolution correction is included. 




Fig. 21. Apparent flux for the same objects in different 
IR/submm wavebands, as a function of redshift for the 
Einstein-de Sitter cosmology. Coding for the lines is the same 
as in Fig. n^ 



0.01 Gyr for the bursts, though part of the scatter can be inter- 
preted as larger ages. 

The synthetic SEDs that are produced illustrate the large 
range of IR/optical luminosity ratios. A striking but obvious 
result is that the optical spectrum of a LIRG or ULlRG-type 
galaxy is extremely sensitive to the details of the dust distribu- 
tion. This makes the predictions on the optical counterparts of 
high-z dusty objects difficult to assess. 

We have gathered a small sample of nearby galaxies that 
have been extensively observed at several wavelengths in the 
optical/lR/submm/radio range. The fit of the overall continuum 
with our SEDs provides us with a way of interpolating between 
the data points under assumptions that have a physical mean- 
ing. We then generate a sequence of galaxies with IR luminosi- 
ties that span several orders of magnitude. The spectral trend is 
very similar to what has been shown by e.g. Sanders & Mirabel 
(1990 on a smaller sample. We have established that there is 



Fig. 20. Same as Fig. n9^ for a flat cosmology where Ho = 65 
km s^i Mpc^i, Oo = 0.3 and Oa = 0.7. 



a very good correlation between the star formation rate and the 
bolometric infrared luminosity, and we have calibrated the re- 
lation. In contrast, the UV fluxes in these objects are strongly 
affected by extinction. Template spectra of the 17 objects listed 
in table H^ generated with STARDUST, are available upon e- 
mail request ( devriendgiap . f r or guider@iap . f r ) 
and can be used for studies of local and distant samples. 

We defer the study of high-z sources to forthcoming pa- 
pers. These SEDs have been explicitly designed to be imple- 
mented into SAMs of galaxy formation and evolution where 
the basic free parameters (i,, t, and the gas column densities) 
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Fig. 22. Same than Fig. 21 for the flat cosmology with a non- 
zero cosmological constant. Coding for the Unes is the same as 
in Fig. Ill 



can be computed from general assumptions of the physical pro- 
cesses ruling galaxy formation. The remaining uncertainties are 
due to the possibility that dust properties at high z can scale dif- 
ferently with metallicity, that perhaps the IMF is not constant, 
and that gas inflows/outflow s can strongly affect galaxy evolu - 
tion. In a companion paper ( Devriendt and Guiderdoni, 1999 ), 
we will implement all this in a simple but explicit cosmologi- 
cal framework to try to derive global properties of galaxies and 
trace their evolution. But the final goal is to work with more 
realistic codes of galaxy formation and evolution, especially 
those which take into account the merging histories of dark 
matter halos (e.g. Kauff mann e t al. ( |1993| ), Baughet al. ( |1996| ), 
Somerville and Primack ( |1999| ), Ninin et al. ( |1999| )). 
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